import sys
sys.path.append('/home/aistudio/external-libraries')

import meshio
import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import Delaunay
from meshpy.tet import MeshInfo, build, Options

def convert_and_remove_duplicates(array):
    # Convert array of int32 to list
    array_list = [list(arr) for arr in array]
    # Flatten the 2D list
    flattened_list = [item for sublist in array_list for item in sublist]
    # Remove duplicate elements
    unique_list = list(set(flattened_list))
    return unique_list

def remove_duplicate_points(points):
    unique_points = np.unique(points, axis=0)
    return unique_points

def remove_duplicate_rows(arr):
    sorted_rows = [tuple(sorted(row)) for row in arr]
    unique_rows = np.unique(sorted_rows, axis=0)
    return np.array(unique_rows)

def plot_points_xy(points, color="k", marker="."):
    figure = plt.figure()
    axis = figure.add_subplot(111)
    axis.scatter(points[:, 0], points[:, 1], s=0.1, color=color, marker=marker)
    plt.axis("scaled")  # 设置x轴和y轴相同的缩放比例
    
def plot_points_3d(points):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(points[:, 0], points[:, 1], points[:, 2], s=0.1)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    # 显示图形
    plt.axis("scaled") 
    plt.show()

def find_left_points(arr1, *arrs):
    # 创建布尔索引
    mask = np.ones(len(arr1), dtype=bool)
    for i, row in enumerate(arr1):
        for compare_row in np.concatenate(arrs):
            if np.array_equal(row, compare_row):
                mask[i] = False
                break
    # 使用布尔索引筛选出第一个数组中没有出现的元素
    result = arr1[mask]
    return result

def find_points_on_truck(all_mesh_points):
    x_min, x_max = 0, 4
    y_min, y_max = 80, 90
    z_min, z_max = 18.5, 20   
    boundary_points = []
    for point in all_mesh_points:
        x, y, z = point
        # 检查是否在长方体的边界上
        if np.isclose(x, x_min) or np.isclose(x, x_max):
            if y_min <= y <= y_max and z_min <= z <= z_max:
                boundary_points.append(point)
        elif np.isclose(y, y_min) or np.isclose(y, y_max):
            if x_min <= x <= x_max and z_min <= z <= z_max:
                boundary_points.append(point)
        elif np.isclose(z, z_min) or np.isclose(z, z_max):
            if x_min <= x <= x_max and y_min <= y <= y_max:
                boundary_points.append(point)
    return np.array(boundary_points)

def convert_small_numbers_to_zero(arr, threshold = 1e-8):
    condition = np.abs(arr) < threshold
    arr[condition] = 0
    return arr

'''----------------------------- Read the .msh file -----------------------------------'''
mesh = meshio.read('original_half_truck.msh', file_format='ansys')

# # Get the node coordinates and cell information
points0 = convert_small_numbers_to_zero(mesh.points)
print("从fluent原网格中读取节点数：", points0.shape[0])

all_points = remove_duplicate_points(points0)
print("删除重复点后剩余节点数：", all_points.shape[0])

# # 内边界点、外边界点和内部点(内外边界点有重合的点)
inner_points = find_points_on_truck(all_points)
print("汽车表面的网格节点数：", inner_points.shape[0])
# 删除外边界中重合的点
other_points = find_left_points(all_points, inner_points)
# 拼接
all_mesh_points = np.vstack((inner_points, other_points))

# 四面体化，生成内外两个凸包，然后合并并删除重复的面
outer_tet = Delaunay(all_mesh_points, furthest_site=False, incremental=False, qhull_options=None)
inner_tet = Delaunay(inner_points, furthest_site=False, incremental=False, qhull_options=None)

wrong_convex_hull = np.vstack((inner_tet.convex_hull, outer_tet.convex_hull))
right_convex_hull = remove_duplicate_rows(wrong_convex_hull)

mesh_info = MeshInfo()
mesh_info.set_points(outer_tet.points)
mesh_info.set_facets(right_convex_hull)
mesh_info.set_holes([(2, 85, 19.25)])
mesh = build(mesh_info, options=Options("p"))

# 新网格的信息
nodes = []
for i, n in enumerate(mesh.points):
    nodes.append(n)
nodes = np.array(nodes, dtype=int)

cells = []
for j, c in enumerate(mesh.elements):
    cells.append(c)
cells = np.array(cells)

print("\nDelaunay四面体化完成!")
print("\nNew mesh infos:")
print("-Nodes:", nodes.shape[0])
print("-Cells:", cells.shape[0])

# 保存vtk文件
vtkFileName = "truck_tet_remesh.vtk"
mesh.write_vtk(vtkFileName)

# 使用meshio读取VTK文件
# mesh_vtk = meshio.read(vtkFileName)

# # 将网格对象保存为Gmsh支持的格式（例如.msh文件）
# mshFileName = "truck_tet_remesh.msh"
# meshio.write(mshFileName, mesh_vtk, file_format='gmsh22', binary=False)
print("\n新网格已成功保存为vtk文件:", vtkFileName)
